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& , Abstract. We present a detailed study of the spectral properties of a locally 

S-h ' correlated site embedded in a BCS superconducting medium. To this end the Anderson 

■ impurity model with superconducting bath is analysed by numerical renormalisation 

group (NRG) calculations. We calculate one and two-particle dynamic response 
function to elucidate the spectral excitation and the nature of the ground state for 
different parameter regimes with and without particle-hole symmetry. The position 
and weight of the Andreev bound states is given for all relevant parameters. We 
also present phase diagrams for the different ground state parameter regimes. This 
, work is also relevant for dynamical mean field theory extensions with superconducting 

symmetry breaking. 
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^ . 1. Introduction 

r**- ■ As described by Bardeen, Cooper and Schrieffer (BCS) pQ electrons in condensed 
matter with an attractive interaction assume a superconducting state below a critical 
temperature, referred to as BCS state. In this state electrons with antiparallel spins form 
singlet bound states (S = 0) known as Cooper pairs. This pair formation is a fermionic 
^ ■ many-body phenomenon as it relies on the existence of a Fermi surface [2]. A singlet 
ground state due to many-body effects also occurs in a quite different situation, when a 
magnetic impurity is embedded in a metallic host [31 H] . This state, known as a Kondo 
singlet, occurs because the electrons in the metal at low temperature experience a large 
effective coupling to the localised impurity spin. As a consequence it is energetically 
favourable to screen the local moment, resulting in a (Kondo) singlet state (S = 0). 

The BCS superconductivity and the Kondo effect, are important topics in their own 
right, and have been extensively studied by the condensed matter physics community. 
The interplay and competition of these two effects have also attracted a lot of interest 
because metals with magnetic impurities can be superconducting at low temperatures 
EJ [7J [8j [10] • The problem of dealing with the two effects together is complicated 
because the magnetic impurities have a disruptive effect on the BCS superconducting 
state and the Kondo singlet formation leads to a breaking of the Cooper pairs. For 
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a recent review on this topic we refer to [11] and references therein. Here we address 
a particular aspect of the problem which has not so far received much attention, the 
effects of the superconductivity on the local spectral properties of the impurity. As in 
earlier studies, we take the BCS superconductor as a fixed reference system and take 
as a model for the impurity an interacting Anderson model. We employ the numerical 
renormalisation group method (NRG), which is a reliable approach to calculate low 
temperature spectral functions. 

From earlier studies of this model, we know that if the interaction U at the impurity 
site is weak, the ground state is dominated by the superconducting behaviour and the 
singlet is predominantly a superconducting one. However, if there is a strong repulsion 
at the impurity site, such that single occupation is favoured, we have a situation where 
a single spin is coupled to the superconducting medium. If the superconducting gap A sc 
is very small then, similar to the case with a normal, metallic bath, the ground state 
is a singlet, more specifically a Kondo singlet. If this gap is increased, however, it is 
not possible to form a Kondo singlet, due to the lack of states in the vicinity of the 
Fermi level, and the ground state becomes a doublet (S = 1/2), corresponding to an 
unscreened spin at the impurity site. This ground state transition at zero temperature is 
an example of a quantum phase transition which occurs for a level crossing that depends 
on a system parameter [12J. The relevant energy scales for this singlet-doublet transition 
to occur in the Kondo regime are the Kondo temperature Tk and the superconducting 
gap A sc . There have been numerical renormalisation group (NRG) studies for the Kondo 
model [131 E] an d Anderson model [15] with superconducting bath. In these works the 
estimate for the ground state transition is given by T K /A SC ~ 0.3, i.e. for T K /A SC > 0.3 
we have a singlet ground state (S = 0) whilst for Tk/A sc < 0.3 the ground state is 
a doublet. We can also consider the transition for a fixed value of A sc and increasing 
values of the local interaction U . In this U increases in the local moment regime, 

Tk decreases until the singlet to doublet transition occurs at a critical value U = U c . 

Due to the proximity effect there is an induced symmetry breaking on the impurity 
site. As a consequence localised excited states (LES) inside the superconducting gap 
can be induced at the impurity site. Such states are well known from superconductor- 
normal-superconductor (SNS) junctions and are usually called Andreev bound states. 
For a weak on-site interaction the ground state of the system is usually a superconducting 
singlet (S = 0) and the LES is an S = 1/2 excitation. It is found that at the ground 
state transition the bound state energy of the LES becomes zero as measured from the 
centre of the gap. This is related to the fact that the level crossing occurs there. 

In recent years detailed measurements on quantum dot structure have enabled one 
to probe strong correlation effects [161 EJ- In these experiments a quantum dot is 
coupled to two leads, which can be superconducting. In such situations finite voltage 
induced currents p2J [19j EHl El] and Josephson currents [22], induced by a phase 
difference, were observed experimentally. For a theoretical description of this situation 
it is important to characterise the Andreev bound states in the gap accurately. Many 
of the more recent theoretical work [231 1211 [251 [261 I2H1 [29], focus on a quantum dot 
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embedded in two superconducting baths with different (complex) superconducting order 
parameters. These situations with two channels and with Josephson or nonequilibrium 
currents will, however, not be covered in this paper. 

For the analysis presented here, which focuses on the spectral properties of locally 
correlated electrons in the superconducting bath, we use the NRG approach. We start 
by outlining some of the details of the NRG calculation with a superconducting medium 
in section 2. We also describe an analysis of the Andreev bound states in the gap in 
terms of renormalised parameters, and discuss the limit of a large gap. In section 3 we 
present results first for the model with particle-hole symmetry. For low energies within 
the superconducting gap we calculate the position and weight of the LES and also give 
the values for the induced anomalous on-site correlation. We also present singlet-doublet 
ground state phase diagrams for the symmetric and non-symmetric cases. The study 
is based on numerical renormalisation group (NRG) calculations, which are capable of 
describing the full parameter range from weak to strong coupling reliably. There have 
been a number of NRG studies of this situation in the past p31 El E2 ETj. However, 
the dynamic response function have not been addressed in a satisfactory way. Here we 
present a thorough study of ground state and spectral properties, which will also be of 
interest for cases where the AIM is used as an effective model for superconductivity in 
the dynamical mean field theory (DMFT) framework. 



2. The Anderson model with superconducting medium 

In the following we consider the Anderson impurity model (AIM) in the form 

H = H d + H mix + H sc . (1) 
The local part H d , which describes an impurity or quantum dot, is given as usual by 



(2) 



with the impurity level e d and an on-site interaction with strength U . Also the mixing 
term has the usual form, 

^ mix = E^( c L c ^ + h - c -)- (3) 

We define T = irV 2 p c as the energy scale for hybridisation, where p c = 1/2D is the 
constant band density of states of a flat band without superconducting symmetry 
breaking. The superconducting medium is given in a BCS mean field form 

H sc = £ £fe4 )CT c fe)CT - A sc [ c fc,T c -fc,| + h - c -l ' ( 4 ) 

fc,CT fc 

where A sc is the isotropic superconducting gap parameter, which is taken to be real for 
simplicity. In equation (jlj) the summations runs over all in a wide band. Another 
energy scale cj d , the Debye cutoff in BCS theory, could enter at this stage to restrict 
the summation. As shown in reference [13] with a scaling argument, this effect does 
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not alter the results substantially and merely leads to slightly different parameters. The 
choice here corresponds to uj^ = D, which was also assumed in earlier work [T3] |15]. In 
appendix A we derive the equation for the non-interacting local d-site Green's function 
matrix of the system (1A.9j) . 

2.1. The numerical renormalisation group (NRG) approach 

For the NRG approach we have to derive a discrete form of the Hamiltonian, which 
can be diagonalised conveniently in a renormalisation group scheme descending to lower 
energies. This is done in an analogous fashion as for a metallic medium described in 
[SniET]. Essentially, there are three steps which only affect H m[x and H sc : 
(1) Mapping to a one-dimensional problem, (2) logarithmic discretisation and (3) basis 
transformation. We obtain 



/ 2T 

w^ = y^E(/oU^+h-c.), (5) 

and 

N a N 

H?JD = £ ln+ i(fUn + i,<r + h.c.) - Y.Ulfli + h-c.) (6) 

cr,n=0 n=0 

where the parameters 7„ have the usual form [4j|. For more details we refer to earlier 
work plCEH]. 

The iterative diagonalisation scheme is set up in the same way as in the standard 
NRG case. Due to the anomalous term in the superconducting bath the charge Q is 
not a good quantum number of the system. Thus eigenstates can only be characterised 
in terms of the spin quantum number S. The coefficients 7„ fall off with n, but the 
second term in (jBJ) does not. Thus the superconducting gap becomes a dominating 
energy scale for large n and a relevant perturbation. It does not make sense to continue 
NRG iterations down to energies much below this scale as there are no continuum 
states anymore in the gap. Therefore, we stop the NRG procedure at an iteration 
N = N max , such that the typical energy scale A - ^™* -1 " 2 is not much smaller than 
the superconducting gap A sc . In practice the number of NRG iterations N is between 
20-50 depending on the magnitude of the gap, where we chose A = 1.8 in all cases. We 
usually keep 800 states and the A\ factor [31] is taken into account in the calculations. 

The NRG approach constitutes a reliable non-perturbative scheme to calculate 
T = ground state properties of a local interacting many-body problem. By putting 
together information obtained from different iterations dynamic response functions can 
also be obtained [I]. Here we calculate these spectral functions in the approach [32l [33] 
based on the complete Anders Schiller basis [31]. The Green's function of the interacting 



problem is given by the Dyson equation (jA.llj) . which involves the self-energy matrix 



E(a>). In appendix B we describe how the diagonal part of the self-energy E(u;) = Su(a)) 
and the offdiagonal part of the self-energy H oS (uj) = Y,2i{^) can be calculated from 
dynamic response functions in the NRG calculation, which is in analogy to the method 
described in reference |35|. 
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2.2. The Andreev bound states 



The denominator of the d-site Green's function, equation ( 1A.11I) . can vanish inside the 



gap | a; | < A sc . As the imaginary part of the self-energy is zero in the gap this leads to 
excitations with infinite lifetime there. They correspond to the localised excited states 
(LES) or Andreev bound states. For the non-interacting case they are determined by 
the equation D(u) = [cf. eq. flA.lOD ]. 

c 2 -^-r 2 + ^ = o, (7) 

where the function E(u) is given in equation (jA.70 . The terms in equation ([7j) are 
functions of uj 2 , so if E® is a solution so is — E®. In general, in the interacting case we 
have to analyse the equation 



E M 



TA 



u-e d +-^-E(uj) uj+e d +-^-+Z(-ujy - ^pl-E oa (-aj)* =0.(8) 



E(ijj) 



rA 



EiUj) 



E{u) 

Once the self-energies are calculated it is possible to solve this equation iteratively. 
Here, we will develop a simplified description by using a low energy expansion of the 
self-energy. First note that in the gap, \u\ < A sc , Im£(co>) = Im£ ofT (a;) = 0. We expand 
the real part of the diagonal self-energy £(a>) to first order around uj = 0, which is 
motivated by the Fermi liquid expansions for the normal metallic case and justified by 
the numerical results for the behaviour for low frequency. The offdiagonal self-energy 
is approximated simply by the real constant S (0). This approximation for the self- 
energy is easy to justify if the gap is small parameter, such that it only covers small 
values of uj. The main objective is to present a simplified picture for the analysis of 
the Andreev bound state for the interacting system. We do not expect to be able to 
describe the system near the quantum phase transition accurately like this, and other 
limitations will be seen in the results later. Hence, we find instead of (jSJ) the simpler 
equation 



2T[uj 2 + A sc zE off 



uj 



.2 ~2 t-i2 „2vioff 



r 2 - z 2 z ott (o) 2 + 1 '7 C ; v 7J =o, (9) 



E(uj 



where renormalised parameters e d = z\e d + E(0)] and f = zT were introduced. As 
usual z~ x = 1 — S'(0). Renormalised parameters for the analysis of the Andreev bound 
states were also considered in reference [251 136]. The definition here corresponds to the 
renormalised perturbation theory framework for the AIM introduced in [37]. The form 
of the equations (J7J) and ([9]) is very similar and both can be easily solved numerically 
to give the bound state solutions to = E® = aEb, a — ±. Due to the additional 
offdiagonal correlations induced by the self-energy term S ofT (0), a simple interpretation 
of the interacting theory based on using renormalised parameters id, T in equation (J7J 
for the non-interacting theory is, however, not possible. 

Based on the same idea we can give approximate expressions for the weights of the 
bound states w% by expanding the diagonal part of the Green's function around ui = Eg. 
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We can write the retarded Green's function in the gap near the bound states u ~ ±E b 
as 

G{uj) = ^ + < . (10) 

uj — E b + it] uj — E b + ir\ 

Using the above approximation for the self-energy the weights are found to be 

w? = -E(E b ) M- . (11) 

In a more sophisticated approximation one could consider an expansion of the self- 
energies around the bound state energies E b rather than u — 0. Various things can be 
inferred from expression (1111) . First we note that in the particle- hole symmetric case, 
Ed = 0, w b = w b = w b . The weights are proportional to the renormalisation factor 
z. Since z shows a similar behaviour as in the metallic lead case they decrease with 
increasing interaction U according to (TTTT) . One can easily see that for bound state 
energies close to the gap, \E b \ — > A sc , the weights go to zero, w b — ► 0. One finds pj)] 
that for small U/ttT and A sc /T < 1 we have E b — > A sc , and also for large U/txT the 
bound state energy is close to the gap. Therefore the overall behaviour for w b is given 
in such a case by w b — > for small U, then an increase with U to a maximum and a 
decay again for large U [cf. figure [3] later]. At the ground state transition, where E b = 0, 
the weight shows a discontinuity, and from equation ( TTTT) this requires a jump of the 
self-energy as function of U. 



2. 3. The limit of large gap 

In order to obtain some obtain analytical results it is useful is to consider the case 
where the superconducting gap is a large parameter, A sc — > oo [251 [23 [291 EH] • Then 
the problem essentially reduces to a localised model with an anomalous on-site term 
which is of the order of the hybridisation T. We will write it in the form 

#d°° = £&( c k c <,«r - 1) - r[4 >T 4 a + h.c] + ^(5>4* - 1) 2 , (12) 

where £<j = Ed + U/2. Without interaction this Hamiltonian can be diagonalised by 
a Bogoliubov transformation and the excitation energies = + T 2 are found, 
which lie in the gap as T <C A sc as assumed initially. This gives a direct picture of the 
emergence of the Andreev bound states for large A sc . 

We can discuss the ground state crossover from the singlet to the doublet state in 
terms of the single site Hamiltonian ( 1T21) . First note that the S = 1/2 (doublet) states, 
|t) and ||), are eigenstates of (|T2|) with zero energy. The S = singlet states, empty 
site |0) and doubly occupied site are not eigenstates of (|T2l) . However, the linear 
combinations in the "BCS-form" , 

I* i> = u d \0) + v d \U), |* 2 > = v d \0) - u d \U), (13) 
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are eigenstates with eigenvalues E\ = —E d + U/2 and E 2 = E d + U/2, respectively. The 
coefficients u^, v d are given by 

«S = i(i + |). ^~w} (") 

The ground-state is therefore a singlet as long as E\ < and a doublet otherwise. The 
condition E\ = or 

t2 p2 1 

W + lP = i (15) 

defines therefore the phase boundary for the transition. It is a semicircle in the (£d/U)- 
(T/[/)-plane with radius 1/2, which is shown in figure [12] later. How this phase boundary 
looks like for finite gap A sc will be investigated in section 13.21 when we look at the 
situation away from particle-hole symmetry. In the case of particle-hole symmetry 
£d = and condition (jl~5]) reduces to V — U/2. 

Having established the formalism and the most important relations we will in 
the next section present results for spectral behaviour of the symmetric AIM with 
superconducting bath with a finite gap parameter. 



3. Results 



In this section we present results for the local spectral properties. The diagonal and 
offdiagonal Green's functions are calculated within the NRG framework usually from 
the Lehmann representation, 

PM = \ EK m l c dN| 2 <^ - (E m - E n )](e-^ + e-^), (16) 

m,n 

and similar for the offdiagonal Green's function. As in this procedure the discrete 
excitations for the spectral peaks in the Green's functions have to be broadened, it is 
not straight forward like this to obtain the sharp spectral gap at \u\ = A sc expected for 
T = . As detailed in appendix B, we can, however, determine the self-energy matrix 
from the one-particle Green's function and the higher F-Green's function [cf. eq. ( IB. 41) ]. 
Then we can use the exact expression for the non-interacting Green's function G^(u) 
in equation (1A.9I) . which includes a sharp spectral gap, and the Dyson matrix equation 
(lA.lip to calculate the diagonal and offdiagonal Green's function, G(uj) and G oS (u>) 
respectively. This is the way the Green's functions are calculated for the region outside 
the gap, |u;| > A sc . Inside the gap, \u\ < A sc , we have extracted the weights Wf, and 
positions E° of the delta-function peaks for the Andreev bound states from the NRG 
excitation data for the Green's function directly from the lowest spectral excitation 
(SE) in equation (fi~6j) . These delta-functions are represented by an arrow in the plots. 
Altogether the diagonal spectral function p(u) = —lmG(u))/'K can then be written in 
the form 

P {U) = J2 ^bS(cO - K) + PcontM, (17) 
a=± 
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where p con t(^) is the continuum part for |u;| > A sc . The offdiagonal part of the spectrum 
p oG (uj) = —lmG oS (uj)/7c has a similar general form as the diagonal part, 

p oS (u) = £ - E?) + (18) 

a=± 

where the weights w% can have positive and negative values. For half filling the spectrum 
p {uj) is an asymmetric function of u. 



3.1. Symmetric model 

We first focus on the particle-hole symmetric model, Ed = —U/2, where the ratio U/ttT 
and the parameter A sc are the relevant energy scales. 



3.1.1. Spectral functions for small gap In figured] we show the spectral function ffTTj) 
for A sc = 0.005 for the diagonal Green's function at the impurity site for a number 
of different values of U. Here and in the following we take a fixed value for the 
hybridisation, 7rr = 0.2. All quantities can be thought of as being scaled by half 
the band width D — 1. 




In the plot on the left hand side we give the spectrum over the full energy range. When 
the interaction is increased, spectral weight is shifted to higher energies as the atomic 
limit peaks at ±{7/2 develop . We also observe the beginning of the formation of a Kondo 
resonance at low frequencies. For larger U the Kondo resonance becomes narrower, but 
its formation is suppressed in the very low frequency regime because the spectral density 
vanishes in the gap region — A sc < u < A sc . This is not visible on the scale used in 
the left hand panel of figure [U In the right hand panel of figure [1] we give an enlarged 
plot of the gap region, which shows the delta function contributions from the Andreev 
bound states, where the arrows give the position of the bound state E b and their height 
indicates the spectral weight Wb- It can be seen that the position of the bound state 
changes when we increase the interaction. The weight first increases and then decreases 
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as a function of U, which corresponds to the feature which was interpreted earlier using 
equation (ITT]) . It is generally of interest to see how much spectral weight is transfered 
from the continuum to the bound states, and an overview for this is given in the later 
figure M (right) and[9j Note that the largest value of U shown, is greater than the critical 
U c for the singlet-doublet transition (U c /irT ~ 3.2). In the high energy spectrum there 
is no significant change to be seen in the behaviour, however, at low energies we observe 
the crossing of the bound state energies at uj = at U c . 




In figure [2] we show the offdiagonal spectral function (ITS]) for A sc = 0.005 for a number 
of different values of U. In the plot on the left hand side we show the behaviour for 
the continuum part outside the gap. Notice that the frequency range only extends up 
to uj = ±0.1. We can see a peak close to uj = ±A SC , which is suppressed for larger U 
and changes sign towards the singlet- doublet transition. The behaviour of the bound 
state peaks in the offdiagonal spectrum is displayed on the right hand side of the figure. 
We can see similar features as observed before in the diagonal part, i.e. the weight first 
increases with U and then decreases. If we follow the excitations with the weight of the 
same sign we can see, that at the singlet-doublet transition the bound state levels cross 
at uj = 0. 



3.1.2. Bound state behaviour A more detailed analysis of the behaviour of the bound 
state as a function of U/ttT and the gap in the medium A sc is presented in figure [31 On 
the left hand side we plot the bound state energies ±Eb and on the right hand side the 
corresponding weights Wb- 

We can see that in the non-interacting case the bound state energy for the cases with 
small gap (A sc = 0.001, 0.01) is very close to ±A SC and decreases to zero with increasing 
interaction. For a critical value U c the nature of the ground-state changes from a singlet 
(S = 0) to a doublet (S = 1/2) and at this point Eb = 0. For this transition we can 
think of the positive E^ and negative solution E^ for the bound states as crossing 

becomes finite again and 



at uj = 0. When the interaction is increased further, 
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u/tit u/jtr 
Figure 3. Bound state energies Et, (left) and weights Wb (right) for various U/ttT and 
A sc . Both quantities have been scaled by the corresponding value of A sc ; 7rr = 0.2. 



increases with U. The larger the gap A sc the smaller critical value U c for this ground 
state transition becomes. In the case where A sc is of the order of T - as can be seen for 
the case A sc = 0.06 - the bound state energy Ef, lies in the middle of the gap already for 
the non-interacting case, but otherwise shows a similar behaviour as described above. 

On the right hand side of figure [3] the weight Wb of these bound states can be seen. 
We have marked the position U c of the singlet- doublet crossover point by a symbol 
on the x-axis. The two curves for a value of the gap A sc = 0.001 and A sc = 0.01 
have a maximum for some intermediate value of U which is smaller than the critical U c 
for the ground state transition. This behaviour can be understood from the analytic 
behaviour of the explicit equation fTTTj) derived earlier. For the other curve (A sc = 0.06) 
the weight is maximal for the non-interacting case. In all cases the weight becomes very 
small for large U. Note that we plot the weight scaled by the gap parameter, -u; b /A sc , 
and therefore the absolute values are larger for the cases with larger superconducting 
gap. At the singlet-doublet transition we can see discontinuous behaviour as the weight 
changes sharply. This is a feature of the zero temperature calculation, where the matrix 
elements in the Lehmann sum ( fl6l) change their values discontinuously when the levels 
cross on increasing U, such that the nature of the ground state changes. It can be seen 
for the anomalous correlations (d^di) in figure El later, as well. For finite temperature 
this discontinuity becomes smooth. 

3.1.3. Spectral functions for larger gap In figure H] we show for comparison the diagonal 
spectral function for a larger gap A sc = 0.02 for the diagonal Green's function at the 
impurity site for a number of different values of U. 

The overall picture on the left is similar to the case in figured] with the smaller gap. Due 
to the larger gap the formation of the central Kondo resonance is completely suppressed, 
but the high energy spectrum is as before. From the behaviour within the gap (right 
side in figure @J we can see that the bound state position goes to zero for smaller 
U values than in the case A sc = 0.005, and hence the ground state transition occurs for 
smaller U c for the larger gap (U c /nT ~ 2.03). This was analysed in detail in figure [3] 
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above. For the values of U shown the spectral weight of the bound states w& decreases 
with increasing U. The weight to& of the peaks in the gap has been scaled differently in 
figures [1] and HI so that their height should not be compared directly. 

The spectral function of the offdiagonal Green's function at the impurity site ( fl8l) 
for this value of the gap, A sc = 0.02, is shown in figure [5]for a number of different values 
of U. 




°' 8 -0.4 -0.2 0.2 0.4 0.6 -3 -2 -1 ,0 1 2 3 

Figure 5. The spectral density p oS (cu) for various values of U for the whole energy 
regime (left) and the region in the gap (right); A sc = 0.02 and 7iT = 0.2. 

For larger frequencies outside of the gap (left) we can see a peak near uj = A sc , whose 
height is reduced on increasing U . At larger frequencies we find that the tails develop 
a broad peak for larger values of U. This has not been observed in the case with the 
smaller gap shown in figure EJ Also a sign change of the low energy peak is found as 
before. The behaviour near and in the gap (right) can be understood as before, where 
in this case we have shown two values of U with a singlet ground state and two with a 
doublet ground state. 

3.1.4- Analysis of bound states with renormalised parameters In section [2] we have 
discussed how the bound state energy, which so far was deduced from the spectral 
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excitations (SE), can also be calculated from the bound state equation (BE) (jSJ). 
The latter was derived by expanding the self-energy to first order. It involves the 
renormalised parameters id, T and the constant value of the offdiagonal self-energy 
E off (0). In figure[6]we compare the bound state energies calculated by these two methods 
for two values of the gap A sc = 0.005 (left) and A sc = 0.06 (right). 




u/itr u/jtr 
Figure 6. Bound state energies Ef, as calculated from the spectral excitations (SE) and 
from the bound state equation (BE) ((SJ with renormalised parameters for A sc = 0.005 
(left) and for A sc = 0.06 (right) for various U/irT; irT = 0.2 is fixed. 

We can see that for values of U < U c the agreement is excellent in both cases. However, 
when U > U c we find less accurate values with the method based on bound state equation 
(BE) with renormalised parameters. Since the method to calculate the bound state 
energy from the NRG spectral excitations (SE) is very accurate we expect inaccuracies 
to be found in the BE method. Indeed, the closer inspection of the numerical results 
for the diagonal and off-diagonal self-energies reveals that the linear and constant 
approximation made in section [2721 to derive the bound state equation with renormalised 
parameters (jSj) becomes less applicable for U >U C . The self-energy displays additional 
features there. 

In section [2] we have also derived an expression ffTTl) for the weights Wf, of the bound 
states in the gap. It can be expressed in terms of the renormalised parameters s~d, T, the 
offdiagonal self-energy E off (0) and the bound states energy E b . In figure [7| we compare 
the weights calculated from the spectral excitations (SE) with the ones from the bound 
state equation (BE) analysis with renormalised parameters. We show the results for the 
same parameters A sc = 0.005 (left) and A sc = 0.06 (right). 

We can see for both cases that the overall behaviour of the weights as a function of 
U is described reasonably well by equation (TTTT) . It is, however, clearly visible that 
the agreement is between the SE and BE values is much better in the singlet regime for 
U < U c . This is similar as observed for the values of the bound states energies in figure 
O and the reason for this is the same. The discontinuity for the weight is not reproduced 
by the approximation based on equation ( jTTT) . As can be seen from that equation this 
would require a sudden change in the self-energy as function of U, which was not found 
with sufficient accuracy in the present calculation. This can partly be attributed to the 
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Figure 7. Weights wt, for the Andreev bound states as calculated from the spectral 
excitations (SE) and from the equation with renormalised parameters for A sc = 
0.005 (left) and for A sc = 0.06 (right) for various U/ttT; ttT = 0.2 is fixed. 



broadening procedure involved and to the inaccuracies when calculating the numerical 
derivative. 



3.1.5. Anomalous expectation value and phase diagram The anomalous expectation 
value (d^di) is an indicator for the strength of the proximity effect of the superconducting 
medium at the impurity site and quantifies the induced on-site superconducting 
correlations. In the following figure M we show the dependence of (dj-cfj,) on the 
interaction U/ttT for the same values of A sc as in figure [31 The values are scaled by the 
gap A sc . 




u / % r u/jtr 



Figure 8. Left: Anomalous expectation values (djdj.) as a function of U/ttT for various 
A sc . The values are scaled by the gap A sc ; 71T = 0.2. Right: The total weight of the 
bound states w in the gap as calculated from the spectral excitations as a function of 
U/ttT for various A sc /7rr. 

We see that as a general trend (d^d^) decreases for increasing on-site interaction. This 
is expected since the superconducting correlations are suppressed by the repulsive 
interaction. We have marked the ground state transition with a symbol on the x-axis, 
and we see that (d^di) changes discontinuously in magnitude and sign there. This is 
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characteristic for this zero temperature quantum phase transition. The sign change is 
due to a phase change of n of the local order parameter which occurs at the transition 
as discussed in reference [TT]. In the situation of infinite gap in the medium, which was 
discussed in section |23| (d-frfj.) drops only to zero at the transition point and is zero in 
the doublet ground state. At finite temperature the behaviour becomes continuous. 

An overview of the transfer of spectral weight from the continuum to the bound 
states is shown in figure [8] (right). There we plot the total weight w = w£ + as 
a function of U/ttY for four selected values of A sc /7iT ranging from 0.005 to 1. The 
curves are similar as before in figure [3] and show the discontinuity at the ground state 
transition. Here the values are not scaled by A sc . We can see that the smaller U and 
the larger A sc are the more spectral weight is found in the bound states. In the extreme 
case of A sc — > we have w — 0, and for large gap, A sc — > oo, and small U equation 
( TTll gives w — > 1. The tendency to both of these limiting cases can be inferred from 
figure [8] (right) and we can see that, for instance, for A sc = 7iT already about 80% of 
the spectral weight is in the bound states. 

Summarising the behaviour for different parameters, we present a phase diagram 
for singlet and doublet states for the symmetric model in the following figure [9j 




Figure 9. Phase diagram for singlet and doublet ground-state as a function of A sc /7iT 
and U/irT, where the full line with large dots describes the phase boundary. The dotted 
line corresponds to U/T = 2, which shows the singlet doublet transition for A sc — > oo. 
The dashed line gives the transition as Xk/A sc ~ 0.3 with Tk given in equation (|19p . 
As a background colour we have included the amount of spectral weight transfered to 
the bound states; the discontinuous behaviour at the the singlet doublet ground state 
transitions is slightly blurred in the interpolated representation. 

For small U the ground state is always a singlet. It can become a doublet when U/nY 
is increased. The critical U c for the transition decreases with increasing value of the 
gap A sc as can be seen in the diagram. In the limit A sc — > oo, the critical interaction 
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is given by U c /tiT = 2/ir, which is shown with a dotted vertical line in the figure. As 
mentioned in the Introduction there have been estimates of the phase boundary for the 
singlet and doublet ground state in the strong coupling regime [131 [15] as ^k/A sc ~ 0.3. 
In this case the Kondo temperature is given as in equation (3.9) in reference |15j . 

T K = 0.182^/^e-^ 8r . (19) 

v nil 

We have added a dashed line representing this result which agrees with the ones 
presented here in the strong coupling regime, but starts to deviate for smaller values of 
U. As a background colour we have included in figure [9] how much spectral weight w 
is transfered to the bound states (The value of w is given by the colour bar on the top 
part of the figure.). As noted before in figure M (right) we can see generally that the 
weight is maximal in the region of large gap and small on-site repulsion U. 

At A sc — > the ground-state is a singlet for any value of U as the Kondo effect 
always leads to a screened impurity spin in a singlet formation. For finite gap the nature 
of the singlet ground state can differ depending on the magnitude of U. It can be a 
singlet corresponding to an s-wave pair like in the wave function given in equation (fl3"]) ; 
which is a superposition of zero and double occupation. This is the natural singlet 
ground state for a BCS superconductor. In the strong coupling regime we can, however, 
also have a screened local spin, i.e. a Kondo singlet. The wave function has a different 
form then and consists rather of a singly occupied impurity state coupled to the spins of 
the medium as many-body state. In our NRG calculations it is not easy to distinguish 
clearly this different nature of the singlet ground states and draw a definite line to 
separate them. We can, however, get an indication for what is favoured from the two 
particle response functions in the spin and in the charge channel. In figure [TU] we show 
the imaginary part of the dynamic charge and spin susceptibility, Xc(u) and Xs(u), for 
A sc = 0.005 and a series of values for the interaction U. 




CO CO 



Figure 10. The imaginary part of the dynamic charge (left) and spin (right) 
susceptibility various values of U; A sc = 0.005 and 71T = 0.2. The scale on both 
axes is the same such that the results can be compared well. 

We can see that the peaks in the charge susceptibility exceed the ones in the spin 
susceptibility for zero and weak interaction indicating the dominance of the symmetry 
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breaking in the charge channel, and a ground state of superconducting singlet nature. 
However, for U/nT > 1 the spin susceptibility develops a large and narrow peak at low 
frequency. This signals the importance of the spin fluctuations and low energy spin 
excitations and indicates a ground states of a screened spin. In contrast the decreasing 
peaks in the charge susceptibility for large U is consistent with the effect of suppression 
of the on-site superconducting correlations. 



3.2. Away from particle-hole symmetry 

So far we have considered the special situation of particle-hole symmetry, = —U/2. 
In this section we will briefly discuss a few aspects that change in the situation away 
from particle- hole symmetry. Let us consider the case where for a given gap A sc , on-site 
interaction U, and hybridisation T, the ground-state of the system is a doublet at half 
filling, £d = 0. When ^ is increased, we find that a transition to a singlet state can 
occur at a certain value This is illustrated in the following figure [HI where we have 
plotted the bound state energy Ej, for fixed A sc = 0.01, two values of U/ttT = 3, 5 and 
a series of values of the on-site energy scaled by U, £d/U. As before we have 7iT = 0.2. 




The critical interaction for the ground state transition for this case at half filling is 
U c /nT ~ 2.6, such that both cases possess a doublet ground state for ^ = 0. We can 
see that with increasing asymmetry &j the bound state energy \Eb\ first decreases towards 
zero and then increases again in the singlet regime for £<j > As in the symmetric case 
the singlet-doublet transition is accompanied by l-E^I = 0. The weights w b for these 
bound states are shown on the right hand side of figure [TTJ Away from particle-hole 
symmetry the weight w b for the positive energy E b and w b the one for the negative 
bound state E b ~ are not equal, as was already pointed out below equation fTTTj) . We 
can see that the weights wf start to assume different values when ^ is increased from 
0. At the ground state transition the values change discontinuously similar as observed 
in the half filled case. If we follow both the positive weight w b and the negative w b 
separately the weights cross at the transition point. If, however, we think of the bound 
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states as crossing at zero, i.e. «-> at the transition, a more direct connection can 
be deduced from the results shown. In the singlet phase there is a maximum for both 
the positive and the negative bound state weight, more pronounced for w£. 

Also in the asymmetric case it is possible to calculate the bound state position 
Eb from equation ([9]) and the weights from equation (lllj) employing the renormalised 
parameters. We do not show the plots here, but note that the results resemble figures [6] 
and [7] in the respect that they give good agreement in the singlet regime, but deviations 
for parameters where the ground state is a doublet. 

In the following figure [12] (left) we show the dependence of the anomalous 
expectation value (d^di) on the asymmetry scaled by the interaction £d/U for the same 
value of A sc as in figure [TTJ The values for (d^di) are scaled by the gap A sc . For the 




Figure 12. Left: Anomalous expectation values (d-fdj.) for various U/irT, A sc = 0.01 
and 71T = 0.2. Right: Phase diagram showing the regions for singlet and doublet 
ground state as dependent on Y/U and £,d/U for different values of the gap A sc . The 
full semicircular line corresponds to the phase boundary for A sc = oo as discussed in 
equation 1(15]). 

values of U shown, at half filling the system has a doublet ground state and (d^di) is 
negative. First it does not vary much when ^ is increased, but at the transition to the 
singlet ground state we find, as in the half filled case, a jump to a positive value and 
(d^di) increases to a saturation value on further increasing This value is smaller for 
larger U, similar to what has been found in the symmetric case. 

On the right hand side of figure [T2] we present a global phase diagram of the 
parameter regimes for singlet and doublet ground states for the non-symmetric case. 
This representation in the r/C/-£d/C/-plane is motivated by the result for the phase 
boundary for the case A sc — > oo derived in section [2~3| equation ffT5l) . The semicircle 
corresponding to this case is shown in the figure together with the phase boundaries 
for some finite values of the gap A sc . These are seen to have a similar form, but the 
boundary decreases to smaller values of Y/U with A sc /7rr. Note that the parameters 
on the line on the x-axis, to which the phase boundary contracts in the limit r — > or 
U — > oo, possess a doublet ground state for /U < 1/2. 
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4. Conclusions 

We have discussed and quantitatively analysed the different forms of behaviour that can 
occur for an interacting impurity site in a medium with offdiagonal symmetry breaking in 
the charge channel. This study is motivated by the experimental situations of impurities 
in superconductors and nanoscale quantum dot systems with superconducting leads. In 
the local spectral functions we found that the low energy spectrum is dominated by the 
superconducting gap, and we saw that the lowest excitations in these cases are Andreev 
bound states within the gap region. For higher energies the spectrum resembles the 
form usually found in a metallic bath with broadened atomic limit peaks for large 
U /ttT. The formation of the Kondo resonance, whose width is proportional to Tk, is in 
direct competition with the superconducting spectral gap of magnitude A sc . Therefore, 
depending on the ratio of these parameters a screened Kondo singlet or an unscreened 
local moment is observed. 

The lowest spectral excitations, the Andreev bound states within the gap region, 
change position and weight according to the other parameters. These have been analysed 
in detail in both the symmetric and the asymmetric model. We have given a simple 
interpretation of their position and weight in terms of renormalised parameters. It 
turned out that the assumptions for the definition of these were satisfied better in 
the singlet ground state regime. The reason for this should be subject of further 
investigation. In the quantum dot systems currents have been observed involving 
multiple Andreev processes [TH1 ED] ■ It is expected that a quantitative understanding 
of these currents require accurate information about the weight and position of the 
Andreev bound states, which have been provided here. To study the experimental 
situation in detail and to describe the differential conductance dependent on the local 
bound state behaviour can be subject of a separate publication, where also the details 
of the experimental setup are taken into account more carefully. 

The behaviour of the ground state of the system, which can be a spin singlet 
or a doublet, is summarised in the two phase diagrams in figures [9] and [12l For the 
overlapping parameter ranges our results for the ground state and the locally excited 
states are in agreement with earlier NRG studies [131 HH EH]- Differences can be seen 
in the spectral representation of the bound states in the gap. Here we report delta 
function peaks, whereas an earlier study [27J presented broadened peaks. The method 
of calculating spectral functions and the self-energy used and explained in the appendix 
of this paper will be relevant for extensions of the calculation to the lattice model within 
the dynamical mean field theory framework. There an effective Anderson impurity 
model could be used to study the phases with superconducting symmetry breaking for 
instance in the attractive Hubbard model. 
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Appendix A. Relevant Green's functions 

For the Green's functions it is convenient to work in Nambu space, C\ = (<4-pC d |), 
with 2x2 matrices. The relevant retarded Green's functions are then 

r (,,\ - ttr -rH - ( (( c d,v c h))» (iw^iK \_( G nH \ . \ , , 

' \ c &,V c d,\))u {{Cd,v c d,i))u J \G 21 {lu) G 22 {uj) 
In the NRG approach we calculate Gu and G 2 i directly and infer G 22 (uo) = -Gn(-w)*, 
which follows from G^bH = -G^-w) and Gf^ dv (cu) = -G^|J v (-w)* for 
fermionic operators A, B. Similarly, we can find Gi 2 {us) = G 2 i(—u)*. In the derivation 
one has to be careful and include a sign change for up down spin interchange in the 
corresponding operator combination. 

In the non-interacting case we can deduce the d-site Green's function matrix exactly. 
To do so rewrite the term H sc by introducing the vector of operators and the symmetric 
matrix 

Then H sc can be written as 

H sc = Y;C{A k C k . (A.3) 
k 

The matrix Green's function in the superconducting lead is then given by g k {ioj n ) = 
(iu n l 2 - Ak)- 1 , 

g k {iu n )~ l = iu n l 2 - e k r 3 + A sc ti, (A.4) 

where Tj are Pauli matrices. It follows that 

. iu n l 2 + e k r 3 - A sc ri 

In the wide band limit with a constant density of states the hybridisation term takes 
the form 

v jy^ M = - r E(iu n ) ■ (A - 6) 
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We are mostly interested in the limit of zero temperature here, and the function in the 
denominator E(z) after analytic continuation reads 



e(u) = { -^ MV" 2 - A2 c fOT M > A - . (A . 7) 



'A 2 C - u 2 for \u\ < A ; 
In the non-interacting case for T = 0, we have therefore 



G^r 1 = *1 2 - e d r 3 + r " l2 + ^ cTl . (A.? 



The Green's function is obtained by matrix inversion, which yields 

- - - - rA 



DM 



where the determinant, D(u) := det(G^(uj) x ) is given by 



D(cu) = uj 



1 + 



r 



E(u) 



r 2 A 2 



"-4 (A.10) 



E(lu) 2 

The full Green's function matrix G d (c<j) _1 at the impurity site is given by the Dyson 
matrix equation 

G^r^G^H-EH, (A.ll) 
where we have introduced the self-energy matrix £(u;). 

Appendix B. Self-energy using the higher F-Green's function 

As described by Bulla et al. [35] there is a method to calculate the self-energy employing 
a higher F-Green's function, and it can also be used for the case with superconducting 
bath. In order to derive the equations of motions for the correlation functions, the 
identity 

^((A] B)) u + (([H, A],B)) W = ([A,B} V ) (B.l) 

(rj = + for fermions) is useful. The calculation taking into account all offdiagonal terms 
yields the following matrix equation 

G 1 (u)G d (cu)-UF(cu) = l 2 , (B.2) 

with the matrix of higher Green's functions F_(u;), 

- l ; V F 2l {u) F 22 {uo) J 1 ; 

We have introduced the matrix elements F u (uj) = ((c d ^n^; Fi 2 (uj) = 

((cdA n i,Cd,i))^ F 2i(^) = -((4,i n T;4,T»^ and F 2 2 (^) = -((4,i n r; c d,i))^- In the 

NRG we calculate F\\ and F 2 \ and the others follow from Fi 2 {u) = —F 2 i{—uj)* and 
F 22 {uj) = Fu(-uj)*. We can define the self-energy matrix by 

= UFiu)^)- 1 . (B.4) 
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The properties of the Green's function and the higher F-Green's function lead to the 
relations £ 12 (u;) = S 2 i(— to)* and X^O^O = ~^n(~ uj)* for the self-energies. We can 
therefore calculate the diagonal self-energy E(co>) = Su(w) and the offdiagonal self- 
energy S ofr (c<j) = T,2±(uj) and deduce the other two matrix elements from them. With 
the relation ( 1B.4I) between G, F_ and S the Dyson equation (1A.11I) is recovered from 
(IB.2j) . Therefore, the Green's function can be calculated from the free Green's function 
as given in flA.91) and the self-energy as calculated from flB.41) . This scheme will be 
useful for applications of dynamical mean field theory with superconducting symmetry 
breaking, where the self-energy matrix has to be calculated accurately to find a self- 
consistent solution. 
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